!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Cold bubble test case as in [Straka, Wilhelmson, Wicker, Anderson,
!! Droegemeier, 1993].
!!
!! \n
!!
!! This is the classical test case discussed in <a
!! href="http://dx.doi.org/10.1002/fld.1650170103"><em>J. Straka, R.
!! Wilhelmson, L. Wicker, J. Anderson, and K. Droegemeier, Numerical
!! solutions of a non-linear density current: A benchmark solution and
!! comparisons, Int. J. Numer. Methods Fl., 17 (1993), pp.
!! 1–22</em></a>.
!<----------------------------------------------------------------------
module mod_straka1993_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_straka1993_test_constructor, &
   mod_straka1993_test_destructor,  &
   mod_straka1993_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_neu,        &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "straka1993"
 character(len=100), protected :: test_description(5)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref

 logical, protected :: mod_straka1993_test_initialized = .false.
 ! private members
 real(wp) :: &
   nu,      & ! kinematic viscosity
   theta_s, & ! potential temperature
   xc(3),   & ! perturbation center
   sigma(3),& ! perturbation size
   dtt        ! temperature perturbation
 character(len=*), parameter :: &
   test_input_file_name = 'straka1993.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_straka1993_test'

 ! Input namelist
 namelist /input/ &
   nu, theta_s, xc, sigma, dtt

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_straka1993_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=3+1+10) :: cformat
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_straka1993_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Straka 1993 cold bubble test, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(cformat,'(a,i1,a)') "(a,",d,"(e11.3),a)"
   write(test_description(3),cformat) &
     '  perturbation center "xc"  = ',xc(1:d),';'
   write(test_description(4),cformat) &
     '  perturbation size "sigma" = ',sigma(1:d),';'
   write(test_description(5),'(a,e10.3,a)') &
     '  perturbation amplitude = ',-dtt,';'

   ! No Coriolis force
   phc%omega = 0.0_wp

   t_ref = 300.0_wp

   mod_straka1993_test_initialized = .true.
 end subroutine mod_straka1993_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_straka1993_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_straka1993_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_straka1993_test_initialized = .false.
 end subroutine mod_straka1993_test_destructor

!-----------------------------------------------------------------------

  !> Viscosity
  pure function coeff_visc(x,uu) result(n)
   real(wp), intent(in) :: x(:,:), uu(:,:)
   real(wp) :: n(2,size(x,2))

    n = nu
  end function coeff_visc

!-----------------------------------------------------------------------

 !> Neutral stability atmosphere with temperature perturbation.
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp), parameter :: tpi = 3.1415926535897932384626433832795028_wp
  integer :: l
  real(wp) :: r, dthetal, dttl
  real(wp), dimension(size(x,2)) :: z, theta, pi, p, t

   ! first define a neutral atmosphere: \theta = const
   z = x(size(x,1),:) ! vertical coordinate
   theta = theta_s
   pi = 1.0_wp - phc%gravity/(phc%cp*theta_s) * z

   ! now perturb the potential temperature (but leave pi unchanged)
   do l=1,size(x,2)
     r = sqrt(sum((  (x(:,l)-xc(1:size(x,1)))/sigma(1:size(x,1))  )**2))
     dthetal = 0.0_wp
     if (r.lt.1.0_wp) then
       dttl = - dtt * 0.5_wp*(cos(tpi*r)+1.0_wp)
       dthetal = dttl/pi(l)
     endif
     theta(l) = theta(l) + dthetal
   enddo

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                     ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + phc%gravity*z) ! e
   u0(3:,:) = 0.0_wp                            ! U
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Neutral stability atmosphere
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp), dimension(size(x,2)) :: z, theta, pi, p, t

   ! first define a neutral atmosphere: \theta = const
   z = x(size(x,1),:) ! vertical coordinate
   theta = theta_s
   pi = 1.0_wp - phc%gravity/(phc%cp*theta_s) * z

   ! then translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                     ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + phc%gravity*z) ! e
   u0(3:,:) = 0.0_wp                            ! U
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_outref

!-----------------------------------------------------------------------

 !> The energy flux is deduced by the constant stability atmosphere,
 !! the momentum flux is set to zero.
 !<
 pure function coeff_neu(x,u_r,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: h(size(x,1),1+size(x,1),size(x,2))

  real(wp), dimension(size(x,2)) :: z, pi, rho

   h = 0.0_wp

   ! vertical component of the energy flux
   z = x(size(x,1),:) ! vertical coordinate
   pi = 1.0_wp - phc%gravity/(phc%cp*theta_s) * z
   rho = phc%p_s/(phc%rgas*theta_s) * pi**(1.0_wp/phc%kappa-1.0_wp)
   h(size(x,1),1,:) = nu * phc%gravity * rho

 end function coeff_neu

!-----------------------------------------------------------------------

end module mod_straka1993_test

